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A variant of the Direct Simulation Monte Carlo method is used to study the behavior of a granular 
gas, in two and three dimensions, under varying density, restitution coefficient, and inelasticity 
regimes, for realistic vibrating wall conditions. Our results agree quite well with recent experimental 
$_( , data. The role of the energy injection mechanism is discussed, as well as the behavior of state- 

On' functions, such as pressure, under realistic boundary conditions. Upon a density increase, we find 

\ signals of a clustering transition. 

. PACS numbers: 81.05.Rm, 05.20.Dd, 45.05,+x 

J* ; I. INTRODUCTION 

on ; 

The study of granular systems (GS) is one of the areas of modern physics that has presented the fastest growth lately 
jl| . The interest in such systems ranges from purely theoretical, as is the case of the self-organized critical sand- pile 
model Q , to commercial, as in the applications oriented to the pharmaceutical industry Q . GS are composed of solid 
matter in the form of irregularly shaped macroscopic grains. Their dynamical and statistical properties can also be 
affected by the presence of an interstitial medium, such as air or a liquid. 

Three recent experiments dealt with two-dimensional granular gases Q subject to periodic shaking It was 

observed that certain quadratic degrees of freedom, only indirectly coupled to the driving walls through the phase 
space randomization occurring during collisions, present non-Gaussian distributions in the steady state. Some of 
these distributions were nicely fitted by stretched exponentials ||,8,^| with an exponent a — 3/2, while others deviate 
altogether from the Gaussian behavior in the dense gas regime [7|. Non-Gaussian tails are indeed predicted by the 
kinetic theory of dissipative granular gases ]To| . Thus, these observations indicate that the dynamics of granular 
systems may depart strongly from that of ordinary gases which are well described by Gaussian distributions. 
. In this paper we study how inelasticity and boundary conditions (either energy feeding or absorbing) affect the 
development of the steady state of two and three dimensional granular gases. By employing numerical simulations 
with realistic boundary conditions, we make contact with the experimental settings [^-Q, studying their similarities 
1 with respect to the energy transfer processes. Our numerical work is based on a variant of the Direct Monte Carlo 
Simulation (DSMC) method p!l| . The implementation of this method is simple and yields results in very good 
£^ ■ agreement with the experiments with essentially no fitting parameters. 

We find that velocity distributions and their moments depend strongly on the energy feeding mechanisms, as well as 
on boundary conditions. In our model, energy enters the system through vibrating walls and exits through collisions, 
either among grains, or between grains and fixed walls, which act as energy sinks. For the inelastic granular gas, we 
find that the velocity distributions loose their Gaussian character, as expected. (Recall that only for elastic systems 
the Gaussian form may be associated with the fact that the kinetic energy, a quadratic form of the velocities, is 
conserved.) In the steady state, we also observe that the ratio between mean square velocities along perpendicular 
directions is not unit. Nevertheless, this ratio is found constant over several decades of vibrating wall velocities. As 
' we show below, this is a consequence of energy entering the system from a privileged direction. Finally, for any 
, inelasticity coefficient value, as the density of the gas is increased, an instability towards the formation of clusters 
appears. While our simulations arc not quantitatively reliable in the high-density limit, they do provide a way of 
characterizing the tendency to clusterization, intrinsic to dissipative granular systems. 

Other recent DSMC simulation |l2],[l3| have addressed the problem of a granular gas under the influence of a 
gravitational field. In contrast, we use a different form of DSMC simulation, including some aspects of grain-wall 
and grain-grain collisions that were not taken into account before, and focusing on the effect of energy feeding and 
dissipation. In addition, We also study how the lateral walls affect the internal energy balance. Moreover, our 
implementation of multiple grain-collision has an alternative, stochastic form. 

The remaining of the paper is organized as follows. In Section [nj we describe the granular system we will study. 
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In Section 111, we derive a convenient form for the one-particle density distribution used in the model. In Section IV, 
we define the length and time scales relevant to the problem. In Section |v|, we calculate the pair collisio n pr obability. 



In Section VI, we sketch the main steps for the computational implementation of the model. In Section VII, we show 
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in detail the boundary conditions employed in our version of the DSMC simulation. In Section VIII, we present our 
numerical results. In Section [X, we discuss our results and conclude. 



II. DESCRIPTION OF THE SYSTEM 



The main object of our study is the complex behavior of granular gases, especially that induced by inelasticity and 
fluctuations. The simplified models we are going to analyze are systems composed of N 3> 1 smooth, spherical grains 
of the same diameter d and mass m. We suppose the system to be contained in a two-dimensional rectangular box of 
lengths L x and L y (later we will generalize the model to three dimensions). Rolling, as well as static friction between 
the grains and the two-dimensional plane are neglected. 

For smooth grains, the interaction only happens upon collisions and in the radial direction. This force is composed 
of a conservative, Hertzian, elastic term and a radial friction one fll4| , 

Fi 2 (ri2,ci 2 ) = k (r 12 - d)i i 12 - fi, lc i(ri2, C12 • fi 2 )fi 2 , (1) 

where r 12 is the vector connecting the centers of mass of grains 1 and 2 and c 12 is their relative velocity. The friction 
term is responsible for the loss of mechanical energy that occurs during a granular collision. In general, that loss is 
well characterized by a velocity- and material-dependent coefficient of restitution e(ci 2 • fi 2 ). 

Due to the continuous collisional loss of energy, one needs to inject kinetic energy into the system to prevent it 
from collapsing. In the experiments [0-0, this is done by means of rapidly oscillating walls, with frequency v and 
amplitude A, at opposite extremes of the box. Here we take the limits v — > oo and A — > 0, keeping Av constant. That 
corresponds to the regime exploited in the experiments. In this manner, energy is given to the system incoherently 
through collisions with the moving walls. Thus, the system reaches a steady state dependent only on vA, and e. In 
this case, the granular gas can be characterized by its phase-space density distribution. 



III. DENSITY DISTRIBUTION 



A. Discretization of the distributions 



In general, we do not have the complete information about all positions and velocities of the system. It is usually 
more appropriate to assume that the one-particle distribution characterizes a stochastic process which can be better 
described by M > N "quanta" of density. These quanta will represent a smooth distribution in phase space. The 
smooth one-particle density distribution can be written in terms of the M quanta as 

/(r,c,i)(Ac)V = a(c,Ac;r,a 2 )/ . (2) 

The normalization for the two-dimensional distribution is then given by 



f d 2 r [ d 2 c/(r,c,t) = Y, Z(r,c,i)(Ac)V, 



phase space 

N = M/ , 
* _ N 

* h - w 

where A denotes the area confining the grains and we used the fact that a(c, oo; r, ^4) = M . Thus, a convenient form 
for the two-dimensional /(r, c, t) is 

N 

f{r ' C > t] = Af«W Q(AC ' a) ' (3) 
Similarly, we may write for the three-dimensional gas 

N 

In the following, we proceed to treat the two and three dimensional granular gas in a discretized form. 
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B. Two-dimensional gas 



The total density distribution variation, due to collisions between grains with velocities c' and c, after a time 
interval St is given by 

<5/(c,c') - /(c) /(c')|c' - c|(Ac) 2 d^. (5) 
We may insert Eq. (||) into Eq. (||) and obtain 

Sf(c,c')=a(Ac,a)a(Ac',a) M ^ 4(Ac j 2 ■ (6) 

It is important to remark that the molecular chaos hypothesis for the collisional probability holds well for an inelastic 
granular system in the the gaseous phase . 



C. Three-dimensional gas 

Similarly, for the three-dimensional case we have 

Sf(c, c') = 7r/(c)/(c')|c' - c|(Ac) 3 d 2 St. (7) 
Inserting Eq. (g) into Eq. @ we obtain 

r „ is \ / a / ,N 2 nd 2 \d -c\5t 

Sf(c, c') = a(Ac, a) a(Ac', a) M2 J (Ac)3 ■ (8) 



IV. NATURAL SCALES AND DIMENSIONLESS UNITS 

There are two important length scales in this problem: the grain diameter d and the box length L (for simplicity, 
let us assume a square box hereafter) . We choose to express all lengths in units of d. We also choose a suitably short 
time interval St as our time scale. This St will correspond to the computational time step in real, physical, terms. All 
lengths are scaled by d and all times by St in such a way that (the asterisks denote rescaled quantities) 

St^ D 



f(v,c,t)=r(v*,c*,t*) 1^ 1 , (9) 
where D is the spatial dimension. In terms of dimensionless variables, Eqs. (0) and (H) become 



f*(r*,c*,t*) = — S- „ a(Ac',fl'), (10) 
Ma* 2 {Ac*) 2 K 



and 

7-2 1„/* 



S 2 r(c*,c>*) = a(Ac*,a*)a(Ac>*,a*) M2 ^ {A ^ 2 , ' U I 



respectively. Similarly, the three-dimensional Eqs. (Q) and (||) become 



f*(r*,c*,f) = — ffr „ a(Ac',a'), (12) 
Ma* 3 (Ac*) 3 v ; ^ y 



and 



* 3 /* (c*, c'*) = c*(Ac*, a*) a(Ac>*, a*) , < L*i 
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V. CALCULATION OF COLLISION PROBABILITIES 



Another way of obtaining Eqs. ([ll]) and ( fl3| ) is to calculate the variation of /*(r, c, t) due to the collisional dynamics 
inside phase space regions of volume a 2 (Acp and a 3 (Ac) 3 , respectively. That variation happens in small jumps due 
to the collisions that actually happen between all pairs of quanta. The total change, in two dimensions, is given by 

£ 2 /(c,c') = -a(Ac, a) a(Ac' , a) M J^ c)2 pl%- (14) 

A comparison between Eqs. ( |l4| ) and (](]) yields the pair collision probability 

2D A|c'*-c*| 
Pcoi.st- Mq * 2 (15) 



in dimcnsionless units. Similarly, we obtain, for the three-dimensional case, 

D _ Nn\c'*-c*\ 



Notice that the pair collision probability of density quanta goes as M , decreasing as M — ► oo. However, the collision 
rate for the actual system remains constant, since the decrease in p co \ is compensated by the increase in the number 
of colliding quanta, which goes as M 2 , and a decrease of the normalization factor (which varies as M _1 ). 

In the next Section, we introduce the Direct Simulation Monte Carlo (DSMC) model [O] , which, in the appropriate 
limit of N fixed and M — » oo, reproduces the usual kinetic theory model. 



VI. THE DSMC AND TIME EVOLUTION 



In our model, space is discretized, but velocities are not. For two-dimensional systems, we divide the position space 
into boxes of volume (a 3 in three dimensions), and, inside each box, the collisional dynamics for the density quanta 
happens as in the usual low density limit. A streaming density flux is established among neighboring cells. We allow 
for the presence of external fields, such as gravity, by the inclusion of an external impulse given to each grain at each 
time step. In the following, we describe the DSMC PJ and its evolution. 

Our version of the DSMC consists of M density quanta (DQ), each carrying the contribution Ma i^ Ac yi [or, Ma J^ Ac yi 

in the three-dimensional case], located in cells of area a 2 (volume a 3 in three dimensions). They reproduce exactly 
the smooth density distribution /(r, c, i) in the limit M — > oo and a — > for a fixed N. 
The time evolution is obtained by following the steps shown below: 

• Step 1: Running over all cells, one constructs an ordered list of all DQ inside each cell. For every DQ pair in a 
given cell, one draws a random number in order to check whether a collision happens or not, with probabilities 
calculated by Eqs. ( |l5| ) or (|l6j). If the collision happen, both DQ are taken out of the list, their new velocities 
calculated by drawing an impact parameter b (and an additional angle ip, in the three-dimensional case) corre- 
sponding to the collision, and proceeding to the next pair in the list. Any grain collide no more than once per 
iteration. 

• Step 2: The streaming flux of quanta is calculated for each DQ in an cell. For instance, let us take a quantum 
located at position (k x ,k y ), with k x ,k y integers, and velocity (c x ,Cy). Thus, the distribution function can be 
written as 

/ = / (^X: k y , c x , c y ,t ). 

Its new position in the x-direction is obtained by drawing a random number rj in the interval [0; 1]: it will be 
k x + [c x 5t/a] if r\ < [c x St/a\, where [z] means the integer part of z; it will be k x + [c x St/a] + 1 if rj > [c x St/a\. 
The analog operation is carried out for the flow along y. The presence of an external field g, such as gravity, 
can be taken into account by adding the impulse 

c — > c + g5t, 

to the particles at this step |p7| . 

• Step 3: The walls' boundary conditions are applied (see below). 

After completing step 3 one returns to step 1. In the following, we explain in detail the computer implementation of 
the DSMC model. 
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VII. ALGORITHM FOR THE DSMC 



A. Initial procedures 

We start by setting the initial distribution uniform in space and Gaussian on the velocities: In fi(k x , k y , c x , c y ) oc 
-c 2 . 



B. Iteration at the borders — boundary conditions and auxiliary sites 

To incorporate the effect of shaking into the boundary conditions, we only need to consider grains whose position 
and velocity direction indicate the possibility of an eminent collision with one of the moving walls. For instance, let 
— c y < be the grain velocity along y for a grain located in a cell of y = 1 coordinate at a given instant t. If not 
disturbed, such grain may collide with the moving wall located at y = (the bottom wall) within an interval St. If 
the streaming drawing decides in favor of the collision, the following procedure is adopted. We approximate the wall 
motion by a sawtooth oscillation with velocities ±Wf,. The amplitude A and frequency v are related by Wb = Av, 
where we take the limit A — > and v — > oo. Also, we allow for the grain-wall collision to be inelastic by introducing 
the inelastic coefficient e w . 

There are three possible ways for a grain to collide with a moving wall. First, it can collide frontally, and only once, 
when the wall has a positive velocity +Wb- In this case, the grain final velocity is given by 

c+ = Wb(l + (17) 

Second, it can also collide only once when the wall has a negative velocity — Wb- In this case, we have 

c_ = -W b (l + e w )+c y e w . (18) 

Finally, the grain can collide twice: if the wall velocity is — Wb at the moment of the collision and, afterwards, the 
grain velocity is sufficiently reduced, the wall can hit the grain again. The second collision is frontal. In this case, the 
grain final velocity will be 

C-+ = W b {l + e w ) 2 -c y e 2 w . (19) 

We now proceed to calculate the probabilities p-,p + ,p |_ for each kind of collision. Their expressions depend on 

which range of velocities c y falls. Notice that c y is the absolute value of the granular velocity near the wall. Let us 
define the following useful velocities: 

cib = W b , 

C2b = W b {l + l/e w ), 
c 3b = W b {l + 2/e w ). 

The probabilities follow from a straightforward accounting of the possible wall positions when a grain enters the region 
of size 2A centered at the average position of the oscillating wall. Once we determine the type of collision, the next 
step is to calculate the corresponding wall reflection flow. For each — c y , the fraction of grains colliding with the wall 
is given by 9 — c y 5t/a. For cells neighboring the wall, the grain-wall probability collision is 

9 x corresponding probability x y- 

Thus, for the DSMC model, each grain near the wall is tested by the drawing of two random numbers. The first 
determines whether the grain collides with the wall with probability 9. The second determines which kind of colli- 
sion it will be (c_ i+i |_). Due to the infinite frequency of wall oscillations, the grain-wall collisions arc completely 

uncorrelated. The grain-wall collision probabilities are given in the following. 

1. Single head-on collision (0 < c y < cf b ): p + = 1. 

2. Single head-on collision, or a head-tail collision, followed by a head-on collision (ci& < c y < C2h): p+ = 
l(l + ? ) and P _ + = l(l-^). 
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3. Single head-on collision, or a head-tail collision, followed by a head-on collision, or a single head-tail collision 



(c 26 < Cy < c 3b ): p+ = \ (l + p_ = 1 - ^ (l + £), andp_ + = I 

4. Single head-on collision, or single head-tail collision (c y > cj, b ): p+ = « (1 



"6 



and p_ 



1 



We repeat the procedure above for the top wall, changing signs, site coordinate, and wall velocity accordingly. The 
left and right walls are assumed fixed. 



C. Velocity scale 



The computational time-scale is set by choosing 

W b : t = W b /-^. (20) 

for given experimental values of W b> t and d. For the case of recent experiments on two dimensional granular gas 

the wall velocities were of the order of 1 m/s, and d w x 10~ 3 m. For a computational W b * t of the order of 1, we obtain 

St w lCT 3 s. 



VIII. NUMERICAL RESULTS 



The main motivation for using the DSMC model is to explore, in detail, the behavior of dilute granular gases, such 
as those studied in Refs. Below, we list the main results obtained from the DSMC. All quantities presented are 

rescaled according to procedure of Sec. IV (asterisks are dropped). 

When averaging over time, we used data points obtained sampling the time sequence at every 10 to 20 steps, 
according to the system size, in order to assure statistical independence. In what follows, we make M = N . 



A. Vibrated two-dimensional gas 



The scheme is presented in Fig. It mimics the experimental setting of Rouyer and Menon |6[ . Their results for 
the steady-state grain densities profile can be reproduced quantitatively by our numerical simulations when we chose 
numerical values for grain number, aerial volume fraction, and inelastic coefficients equivalent to the experimental 
ones. For instance, Fig. || shows the vertical (longitudinal) and horizontal (transverse) average grain number per cell 
density profiles for the two-dimensional square box. It is important to remark that the wall rms acceleration in Ref. |j] 
was sufficiently strong to render the gravitational field ineffective. We thus have set g = in these simulations. Notice 
that the grain density along y is lower near the moving walls, reaching a maximum in the middle of the box. The 
opposite behavior happens horizontally, although the density variation is less pronounced. The velocity distributions 
P{c x ) and P{c y ) were evaluated using only data from cells of maximum density along y, namely, k y — 5, 6 for 10 x 10 
cells. 
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FIG. 1. Schematic view of a two-dimensional, vertical granular system. Spatial discretization is implemented by dividing 
the square box containing the spherical grains into equal cells of side a. The double vertical arrows indicate the position of the 
moving walls. 
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cell coordinate 

FIG. 2. Grain density distributions for 10 x 10 cells, a = 5, N = 200, e = 0.5, e w = 0.8, W b = W t = 1, and 3 = 0. 

In Fig. ||, we show the velocity distribution along the direction (vertical) parallel to the moving walls, P(c y ), for a 
total grain number ranging from N = 100 to 300. To facilitate the comparison, the velocities are rescaled by their mean 
square value, a y — (Cy) 1 ^ 2 . These curves agree qualitatively with the experimental results, showing the symmetric 
shoulders characteristic of the energy injection mechanism. In fact, for a given grain number, each distribution is 
formed by the superposition of three identical curves: one centered at c y = and two others centered at c y = ±c™, 
where c y increases with the wall velocity and decreases with e w . Notice that the shoulders do not coincide because the 
rescaling factor <j y varies with the grain number. Another interesting feature is the nearly exponential tails. Similar 
results were obtained by Baldassari et al. in another recent DSMC simulation for the horizontal distributions Il3|. 
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FIG. 3. The vertical (longitudinal) velocity distributions for the same parameters as in Fig. |^. The solid lines are only guides 
for the eyes. 



Figure [|shows the horizontal velocity distribution P(c x ) for the same grain numbers of Fig. ^ and a similar rescaling 
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(a x = (c 2 ,) 1 / 2 ). We observe that, as the number of grains (and consequently the average number of collisions per grain 
per unit of time) is increased, the form of the distribution changes smoothly from nearly Gaussian to, approximately, 
an stretched exponential of the form Ae~^^ Cx l a ^ , with a density-dependent exponent a. For N = 200, the exponent 
a 3/2 fits approximately the tails of the distribution (the normalization condition requires that w 0.797 when 
A = 1 in this case). For N = 300, a is certainly smaller than 3/2. This is illustrated in Fig. ||, where the tails of the 
horizontal velocity distributions are presented in a log-log scale. 
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FIG. 4. The horizontal (transversal) velocity distributions for the same parameters as in Fig. g. The exponent a indicates 
exponential (1), stretched exponential (3/2), and Gaussian (2) curves. 
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FIG. 5. The tails of velocity distributions shown in Fig. 



We repeated the numerical simulations for several values of Wb, W t , and ew We verified that the only important 
aspect about the velocity of the vibrating walls is the injection of energy into the system that they are responsible 
for. The wall velocities and wall-grain restitution coefficient set the steady-state average velocity in the gas, as shown 
in Fig. ^, but do not alter the shape of the (rescaled) velocity distributions. 
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wall velocity (W b = W t ) 

FIG. 6. Grain root mean square velocity dependence on wall velocities. The simulation parameters are the same in Fig. 
Units follow the rescaling introduced in Sec. NVK Notice that a x 7^ a v , but their ratio is constant. 



B. Dependence on inelasticity 

While changes in the grain-wall restitution coefficient only rescale the steady state average square velocity compo- 
nents, we found that the grain-grain inelasticity does set the behavior of the velocities distribution tails for sufficiently 
large and dense systems. That effect can be seen in Fig. ^, where we present the distribution of horizontal velocities 
for different values of e. Notice that the distribution starts as Gaussian when the inelasticity is weak, but evolves 
towards a stretched exponential as e becomes substantially smaller than unit. 
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FIG. 7. The dependence on e of the horizontal velocity distribution. The parameters are the same as in Fig 



The influence of the grain-grain restitution coefficient in the standard deviations of horizontal and vertical velocities 
are shown in Fig. Notice that the ratio a y /cr x approaches 1 as the system becomes nearly elastic. In the opposite 
limit, the asymmetry between longitudinal and transverse average square velocities increases as the system becomes 
more inelastic. The equipartition of average kinetic energy between horizontal and vertical degrees of freedom is 
therefore broken. In this case, it is meaningful to define two granular temperatures for the system, namely, 

T x =a 2 x and T y = a 2 y . (21) 

Pressure is defined on the wall as the average momentum transmitted to it by the colliding grains per unit area and per 
unit time. Note that pressure depends on the inelastic properties of the grains and the wall. We thus also define P x 
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and P y as the average momentum transmitted for the horizontal and vertical walls, respectively. The PTN diagram 
for a fixed system area is shown in Fig. [5| Notice that their interdependence is still linear along a given direction. 
However, dynamical effects due to energy injection mechanism appear related to the higher pressure on the direction 
perpendicular to the moving walls. 
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FIG. 8. Grain root mean square velocity dependence on e (inelasticity coefficient). The simulation parameters are same as 
in Fig. ^, with N = 200. Also shown is the ratio between the rms velocity components. Notice that one does not expect 
<Jx/<Jy — 1 in the elastic limit since energy enters from the wall along x. 



We have also investigated a possible dependence of the P{c y ) and P(c x ) on the functional form of e(ci2). The 
numerical results (not shown) indicate that the influence of the detailed dependence of e on en is weak when e is not 
too small. 
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FIG. 9. The ratio between wall pressure P and temperature T along each direction as a function of grain number (fixed 
volume). The data for different restitution coefficients are shown. Other simulation parameters follow those of Fig. |j| All units 
follow the rescaling introduced in Sec. IV 



The equation of state can be studied locally inside the gas volume ]18[. In the steady-state regime, we found that 
the granular temperature and the density vary inside the box, but their product, proportional to the pressure, is 
constant along a fixed direction. This can be seen in Fig. [To], where we have plotted the ratio between the local 
quantity T x x n near the wall and in the bulk. However, we notice that for large inelasticities the constancy of the 
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product within the gas is no longer valid. In that case, while pressure should still be homogenous along a given 
direaction, the value obtained for the granular temperature is likely to be incorrect. This occurs because we have 
ignored the multiple grain-grain collisions at higher densities that should lower considerably the value of the average 
kinetic energy. Therefore, corrections associated with an increase in density are needed. 

Furthermore, in the present model, we ignore the grain- wall tangential friction. The presence of this type of friction 
in real systems leads to cluster breakdown at the walls, thus reducing density instabilities as those observed in our 
model. It also feeds energy into the degrees of freedom parallel to the vibrating walls during grain-vibrating walls 
collisions, helping to reduce the difference between T x and T y . 
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FIG. 10. The ratio between T x x n (the product of local temperature and density) near the wall and in the bulk. While the 
ratio remains very close to unit for quasi elastic systems, it begins deviating strongly from that rule when the system is close 
to an instability. 



C. Clustering transition 



The model fails when densities become too high. In that regime, the particles tend to form clusters at low temper- 
ature and the collisional probabilities used in the DSCM simulation becomes incorrect. In fact, our method becomes 
unstable for high densities or strong inelasticity Such instabilities can be identified in the present simulation by a 
sudden increase of granular density (sometimes of several orders of magnitude) in certain regions of the box, especially 
along the still walls. Thus, while the method does not allow any quantitative description of the clusterized phase, it 
does provide a way of identifying, within the molecular chaos hypothesis, a lower bound for the onset of grain clusters 

EE- 

The estimate an upper bound for the stability of the realistic granular gas goes as follows. If an excessive density 
increase occurs near a still wall (even at values so large that the excluded volume of the particles would be larger 
than the region containing them) , the center-of-mass of the system will be dislocated towards that wall. This provides 
us with a method to investigate the clustering threshold: If the center-of-mass displacement is larger than a typical 
fluctuative displacement, we check for excess density peaks near the walls. If these are present, the model is assumed 
to fail in that region. We thus plot a density-inelasticity phase diagram by marking the upper boundary of stability 



obtained numerically (see Fig. 11). We observe that for any value of the inelasticity, there is a minimum value of the 



density for which a cluster appears 121 
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FIG. 11. The value of e below which an instability occurs for a given TV. Other parameters follow those of Fig 
dashed line is polynomial (cubic) curve fitted to the data points. For N < 200, no instability was found. 



The 



D. The three-dimensional gas 



The generalization of the two-dimensional simulations to three dimensions can be straightforwardly done in the 
DSMC context. We investigated mainly two regimes (for not too small e), namely, low and moderate densities. To 
facilitate the data interpretation, we adopted hard wall boundary conditions and set gravity equal to zero. The moving 
(equal velocity) walls operated in the x direction only. 

For a small number of grains, N = 5000 in a cubic box of dimensions (5 x d) 3 , the transverse velocities followed 



closely a Gaussian distribution, as in the case of an elastic gas. This can be seen in Fig. 12. The density distribution 
along the direction parallel to the moving walls (y, not shown) is similar to the two-dimensional case, having a 
maximum at zero velocity and two lateral "shoulders" whose positions are related to the vibrating wall velocities, 
signaling that they are caused by the energy injection mechanism. 
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FIG. 12. The distribution of velocity components perpendicular to the moving walls (z direction). System parameters are 
10 x 10 x 10 cells, a = 5, Wb = W t = 1, e = 0.5, e w = 0.8, and g — 0. The distribution along x is essentially identical (not 
shown) 



The moderately dense gas case, with N = 10, 000 in a cubic box of the same volume, presents a much more 
interesting behavior. The velocity distributions strongly depart from Gaussian, or even an stretched exponential 
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behavior, as seen in Fig. |j. 

There seems to be no essential difference between the two- and the three-dimensional granular gases. Our method 
allows us to investigate the moderately dense granular gas limit, using pair collision probabilities that are close to 
actual experimental values. We expect to obtain reliable results as long as the probability of clustering is small. We 
stress that velocities distributions seem to change continuously from Gaussian to power-law behavior, as the effective 
dissipation per grain per unit time is increased. 



IX. CONCLUSIONS 



In this paper, we studied a two-dimensional granular gas via a completely stochastic DSMC model. We used realistic 
boundary conditions for energy feeding and energy absorbing walls. We analyzed the behavior of the system for various 
regimes of granular densities and inelasticity. We observed that the shape of the parallel velocities distribution is 
highly dependent on the energy feeding mechanism. The appearance of distinct granular temperatures for the parallel 
and the transversal directions is one consequence of it. 

Our model reproduces experimental data for granular velocities distributions on a satisfactory way, but fails when- 
ever the density or inelasticity take very large values. It also predicts a different set of granular temperatures for the 
parallel and orthogonal directions to the shaking walls. This is a manifestation of the fact that, for granular inelastic 
systems on a steady-state, the equipartion theorem of equilibrium statistical mechanics does not apply. A recent 
experimental work has observed the same phenomenon for a mix of two different granular systems p2| . 

We also studied some of the limits of such model by looking for the breakdown of the gas phase due to the formation 
of clusters. We were able to estimate an upper bound for the densities necessary for the formation of clusters for a 
given inelasticity value. The type of inelastic collapse associated to the formation of clusters reminds us of experiments 
with vibrated granular systems where the bulk of the system gets clusterized on one side of the volume, while the 
other side remains in the fluidized phase B5(l. 

The equation of state of a granular gas and the behavior of the PTN diagram for moderate densities and inelasticities 
were investigated. In principle, the onset of inelastic collapse could be as well estimated by the behavior of the local 
density-temperature product. 

In summary, in this article, we present a computationally easy tool for the study of the complex behavior of granular 
systems and study some peculiar properties (e.g. equipartition breakdown, non-gaussian behavior, etc.) and some of 
its limitations (e.g. clustering formation). 



X. ACKNOWLEDGEMENTS 



We would like to thank CNPq, FAPERJ, and PRONEX for partial financial support. 



[1] H.M. Jaeger, S.R. Nagel, R.P. Behringer, Rev. Mod. Phys. 68 (1996) p. 1259; Phys. Today 49 (1996) p. 32. 

[2] P. Bak, C. Tang, K. Wiesenfeld, Phys. Rev. Lett. 59 (1987) p. 381. 

[3] T. Shinbrot, F.J. Muzzio, Nature 410 (2001) p. 251; Phys. Today 53 (2000) p. 25 

[4] H.J. Herrmann, S. Luding, R. Cafiero, Physica A 295 (2001) p. 91. 

[5] A. Kudrolli, J. Henry, Phys. Rev. E (Rapid Commun.) 62 (2000) p. 1489. 

[6] F. Rouyer, N. Menon, Phys. Rev. Lett 85 (2000) p. 3676 

[7] D.L. Blair, A. Kudrolli, Phys. Rev. E (Rapid Commun.) 64 (2001) 050301. 

[8] W. Losert, D.G.W. Cooper, J. Delour, A. Kudrolli, J.P. Gollub, Chaos 9 (1999) p. 682. 

[9] J.S. Olafsen, J.S. Urbach, Phys. Rev. Lett. 87 (1998) p. 4369; Phys. Rev. E (Rapid Commun.) 60 (1999) p. 2468. 
[10] T.P.C. van Noije, M.H. Ernst, Granular Matter 1 (1998) p. 57; S.E. Esipov, T. Poschel, J. Stat. Phys. 86 (1997) p. 1385. 
[11] G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Clarendon, Oxford, 1994. 
[12] A. Puglisi, V. Loreto, U.M.B. Marconi, A. Petri, A. Vulpiani, Phys. Rev. Lett. 81 (1998) p. 3848. 
[13] A. Baldassarri, U.M.B. Marconi, A. Puglisi, A. Vulpiani, Phys. Rev. E 64 (2001) 011301. 

[14] N.V. Brilliantov, F. Spahn, J.-M. Hertzsch, T. Poschel, Phys. Rev. E 53 (1986) p. 5382; W.A.M. Morgado, I. Oppenheim, 
ibid. 55 (1997) p. 1940; G. Kuwabara, K. Kono, Jpn. J. Appl. Phys. 26 (1987) p. 1230. 



13 



[15] S. Chapman and T. G. Cowling, The Mathematical Theory of Non-uniform Gases, 3rd. edition, Cambridge Press, Cam- 
bridge, 1990. 

[16] R. Soto and M. Mareschal, Phys. Rev. E 63 (2001) 041303. 

[17] This global change affects the grain distribution mean position in the same way as a uniform gravitational field acting on 

a free particle, i.e., (y) = gt 2 /2. Moreover, the second moment grows linearly with time. 
[18] S. Luding, Phys. Rev. E 63 (2001) 042201. 
[19] J.M. Pasini, P. Cordero, Phys. Rev. E 63 (2001) 041302. 
[20] E.L. Grossman, T. Zhou, E. Ben-Nairn, Phys. Rev. E 55 (1997) p. 4200. 

[21] Large hydrodynamic instabilities appear when the system size is larger than a certain critical value. See, for example J.J. 

Brey, M.J. Ruiz-Montero, Comp ut. Phys. Comm. 1 22 (1999) p. 278, and references therein. 
[22] K. Feitosa, H. Menon, preprint (|cond-mat/011139l|). 



14 



